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We present a method for approximating the solution of the three-dimensional, time-dependent 
Gross-Pitaevskii equation (GPE) for Bose-Einstein condensate systems where the confinement in 
one dimension is much tighter than in the other two. This method employs a hybrid Lagrangian 
variational technique whose trial wave function is the product of a completely unspecified function 
of the coordinates in the plane of weak confinement and a gaussian in the strongly confined direction 
having a time-dependent width and quadratic phase. The hybrid Lagrangian variational method 
produces equations of motion that consist of (1) a two-dimensional, effective GPE whose nonlinear 
coefficient contains the width of the gaussian and (2) an equation of motion for the width that 
depends on the integral of the fourth power of the solution of the 2D effective GPE. We apply 
this method to the dynamics of Bose-Einstein condensates confined in ring-shaped potentials and 
compare the approximate solution to the numerical solution of the full 3D GPE. 



PACS numbers: 03.75.Gg,67.85.Hj,03.67.Dg 



I. INTRODUCTION 

Recent advances in laser-control technology have en- 
abled the laboratory realization of Bose-Einstein con- 
densate (BEC) systems subjected to all-optical poten- 
tials which provide strong confinement in a horizontal 
plane and an arbitrary potential within this plane. These 
potentials can be produced by a combination of a hori- 
zontal light sheet combined with a rapidly moving red- 
or blue-detuned vertical laser that "paints" an arbitrary 
time-averaged optical dipole potential in the horizontal 
plane Horizontal light sheets can also be combined 
with vertically propagating beams in specialized laser 
modes, such as Laguerre-Gauss modes, to produce other 
types of novel potentials 0. In addition to providing 
strong vertical confinement and counteracting the effect 
of gravity, the light sheet provides stabilization against 
dynamic excitations of the condensate 0] as well as ther- 
mal phase fluctuations [|J. 

The ability to create and probe quasi-2D BECs in ar- 
bitray 2D potentials is motivated by several areas of cur- 
rent ultra-cold atom research. For example, condensates 
in toroidal and ring lattices can be studied. Stable states 
of multiple vortices and persistent currents can be created 
and studied by stirring the condensate 0, H, @ ■ There 
are proposals for creating ring lattices and for studying 
non-equilibrium phase transitions within this geometry 
@, HI • Toroidal geometries are well-suited for studying 
topological defects that may appear during a rapid cool- 
ing process that produces a condensate [{J, Ho| . There 
is some indication that stirring within a ring-lattice ge- 
ometry can produce a coherent superposition of states 
with different circulation which can lead to a reduction in 



the threshold of the Mott-Insulator phase transition [11 



These systems also offer an excellent finite sized testbed 
for systems of ultra-cold atoms that mimic condensed- 



matter systems [12J. 

Quasi-2D BECs may also provide a convenient plat- 
form for studying systems of ultra-cold neutral atoms 
that are analogs of electronic materials, devices, and cir- 
cuits [l3[. Such systems are called "atomtronic" because 
strongly interacting Bose gases in a lattice potential are 
analogous to "electronic" systems of electrons moving in 
the periodic lattice potential of a crystalline solid. The 
ability to produce arbitrary potentials in the plane of the 
quasi-2D condensate may enable the controlled study of 
novel atomtronic systems. In particular it may be possi- 
ble to produce circuit-like potentials within the plane. 

The behavior of many of the above-mentioned ultra- 
cold bosonic systems can be described using mean-field 
theory. In this case, the governing equation is the time- 
dependent Gross-Pitaevskii equation (TDGPE) [14[l5j . 
This is a partial differential equation in three space vari- 
ables and one time variable whose solution represents the 
wave function of the single-particle orbital that all of 
the condensate atoms occupy. Experiments conducted 
on these systems typically involve releasing the conden- 
sate for imaging. In this case solution of the 3D TDGPE 
becomes a challenging numerical problem due to the vol- 
ume that must be accounted for in simulating the exper- 
iment. 

In this paper, we present a variational approximation 
to the solution of the TDGPE for these quasi-2D systems 
which produces equations of motion whose numerical so- 
lution can be obtained 100 to 1000 times faster than solv- 
ing the full 3D TDGPE. This approximation is based on 
a variant of the standard Lagrangian Variational Method 
(LVM) (T|| in which some of the variational parameters 
are functions of the space coordinates and time while oth- 
ers are only functions of time. The work presented here 
applies this "hybrid" version of the LVM to a quasi-2D 
systems of bosonic atoms. The hybrid LVM was previ- 
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ously applied to a quasi-lD system where only one di- 
mension was weakly confined compared to the two other 
dimensions flT\ . 

This paper is organized as follows. In Section [TT] we 
describe the LVM and its hybrid form and derive the 
approximate equations of motion. We also derive the 
equations that provide the proper variational stationary 
solution in which a condensate is trapped in a confining 
potential. Section IIIII presents a comparison of the so- 
lution of the hybrid LVM equations of motion with the 
numerical solution of the 3D TDGPE for a BEC con- 
fined in a ring-shaped potential. The parameters for this 
example system were taken from an actual experiment. 
Section HVl presents a summary of the work. 



II. TWO-DIMENSIONAL HYBRID LVM 
EQUATIONS OF MOTION 

The condensate wave function of a BEC that is 
strongly confined in one dimension (z direction) rela- 
tive to the confinement in the other two dimensions {xy 
plane) can often be approximated as the product of a 
function of x and y only with a gaussian function of z 
only. In the mean-field approximation, the actual behav- 
ior of the condensate wave function is governed by the 3D, 
time-dependent, Gross-Pitacvskii equation. However, it 
is possible to find equations of motion from which the 
approximate product wave function can be constructed 
at each moment of time using the Lagrangian Variational 
Method. We briefly describe this method next. 



A. The Lagrangian variational method 

The Lagrangian Variational Method provides approx- 
imate solutions to the 3D TDGPE in the form of equa- 
tions of motion for time-dependent parameters that ap- 
pear in an assumed trial wave function. Thus, in the 
standard LVM, the exact solution of the 3D TDGPE re- 
quiring the solution of a partial differential equation in 
three space and one time variable is traded for the so- 
lution of ordinary differential equations in time for the 
variational parameters of a trial wave function of fixed 
functional form. 

The 3D TDGPE can be written as 

= -^ v2 * + Wr)* + 9N |*| 2 *, (1) 

where M is the mass of a condensate atom, g = Airh 2 a/M 
is the interaction strength of low-energy binary scatter- 
ing events with a being the s-wave scattering length, N 
is the number of atoms in the condensate, and Vt rap (r) 
is the external potential. 

The TDGPE is itself a variational equation-of-motion 



and is derived from the following Lagrangian density: 

£[*] = \ih(^i-^^ t ) + ^j 

rj—x,y,z 

+ V trap (r)*** + igV(**) 2 (*) 2 (2) 

where ^ v = d^/drj and rj = x,y,z,t. The associ- 
ated Euler-Lagrangc equation that produces the TDGPE 
with the above Lagrangian density is given by 




The LVM is an approximation method that produces 
an equation of motion for the n time-dependent vari- 
ational parameters, qk{t),k = 1, ...,n, appearing in a 
given trial wave function = "0triai(r; qi, ■ ■ ■ , q n )- The 
equations of motion for these parameters are obtained 
by inserting the trial wave function into the LVM La- 
grangian density, integrating this over the spatial vari- 
ables: 

£lvm(?i (»,•••, ?«(*)) = J d 3 r£ [Atrial (r,t)] , (4) 
and applying the usual Euler-Lagrangc equations: 

d ( (9Llvm\ <9Llvm n • i 
Jt\-^q-)-^q-^- 3 = 1 >-> n - ^ 

This is the standard Lagrangian Variational Method [HJ . 

The LVM can be regarded as having two limits in terms 
of the chosen trial wave function. The first limit consists 
of choosing a trial wave function where the variational 
"parameter" is ^(r,^. This choice enables the varia- 
tional solution to vary in any possible way. As noted 
above, when Eq. ([3]) is applied to the Lagrangian den- 
sity C[^J] to derive the equation of motion, it turns out 
to be the full TDGPE. In the other limit, the trial wave 
function is chosen to have a fixed functional form of the 
spatial coordinates where the time dependence resides en- 
tirely within a set of variational parameters, {Sj(i)}™ = i, 
so that the Lagrangian depends only on these parame- 
ters, Zzlvm [<7i , • • • j Sn]- The shape of this trial wave func- 
tion can only be varied by changing the values of the 
qj. The equations of motion for the qj(t) are ordinary 
differential equations in time and are obtained from the 
usual Euler-Lagrange equations, Eqs. (|5|). It is also pos- 
sible to choose a "hybrid" trial wave function that plots 
a course midway between these two limits. We describe 
this approach now. 

B. The hybrid LVM 

The "Hybrid Lagrangian Varational Method" (HLVM) 
is an LVM in which the trial wave function consists of a 
completely unspecified function of some of the spatial 
coordinates, <f)(x,y,t), multiplied by a fixed function of 
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the rest of the coordinates that also contains some time- 
dependent variational parameters, qi(t), . . . , q n (t). The 
HLVM is expected to apply to systems where there is 
tight confinement in one or two dimensions. The coor- 
dinates appearing in <f> are those for which the confine- 
ment is weak while the trial wave function is assumed to 
be gaussian in the coordinates of tight confinmcnt. The 
HLVM for tight confinement in two dimensions has been 
studied earlier [ItJ ■ Here we study the case where there 
is tight confinement in one dimension only. 

Coupled equations of motion can be derived from a 
"hybrid" Lagrangian which is constructed by integrat- 
ing the Lagrangian density L in Eq. ([2]) over the space 
coordinate(s) of the tightly confined direction(s). The re- 
sulting hybrid Lagrangian can be used to derive coupled 
equations of motion for both <p and the set of variational 
parameters {q\, . . . ,q n }. 

Before proceeding with the derivation of these equa- 
tions of motion, we will first introduce scaled variables 
and rewrite the LVM equations in terms of these vari- 
ables. Scaled units are referenced to a chosen unit of 
length, denoted by Lq, and scaled spatial coordinates are 
given by 



x - V x 
x=—, y=-j~, z=—. 

L,Q LlQ Ijq 



(G) 



Energy and time units are defined in terms of Lq enabling 
the definition of a scaled time: 



h h 
_c/fi o , 1 n = ■ 

u 2MLq Eq 



t 



(7) 



Hereafter barred symbols will denote quantities expressed 
in their appropriate scaled units. It will also be conve- 
nient to express the solution of the 3D TDGPE in terms 
of scaled units: 



V(r,t) = L 3/2 $(r,t). 

In terms of these variables the TDGPE becomes 
,9$ _ /«9 2 $ <9 2 $ <9 2 $ 



(8) 



8t- ' 0f + ^) +V ^ + ~ 9Nm *' 

(9) 

where g = g^E^L^). In scaled units, the Lagrangian 
density takes the form: 



+ i4 r ap$*$ + i<?iV($*) 2 $ 2 , 
and the scaled Euler-Lagrangc equation becomes: 



V — ( — \ - — 

^ dri I d<f>* J ~ d<f>* 







(10) 



(11) 



Now we turn to the description of the hybrid Lagrangian 
Variational Method. 

In deriving the HLVM equations of motion we will as- 
sume that the trapping potential can be written (at least 



approximately) as the sum of a part that depends only 
on the loosely confined coordinates (here x and y) and 
a part that is harmonic in the tightly bound direction. 
Under this assumption we can write the potential as: 



V tra p(x, y, z) = V\\ {x, y) + A 2 z 2 



(12) 



This form of the potential applies in many realistic exper- 
imental cases such as the painted potentials mentioned 
earlier. 

The trial wave function for the HLVM equations of 
motion is written as follows: 

*tn.i(S, = (f>(x,y,t)A(t)e-* 2 / 2 ™ 2 W+0(m\ (13) 

Here the trial wave function is a product of a completely 
unspecified function 4>(x,y,t) with a gaussian function 
having a time-dependent width w(t) and quadratic phase 
coefficient f3(i). These are the variational parameters 
that will appear in the HLVM equations of motion. 

The first step in the hybrid LVM consists of construct- 
ing a "hybrid" Lagrangian by integrating only over the 
spatial coordinate along which the system is strongly con- 
fined: 



^hybrid [<?!>, A, W, 0] 



dzC 



<t>Ae~ 



(14) 



The resulting hybrid Lagrangian is given by 

£hybri<#, W,0\ = § ((j)(/>f - (j)*4>t) + <t>l<t>x + <Pl^y 



+ iA 2 ^ 2 j + igV(r) 2 (0) 2 (^-). 

(15) 

In the above we have eliminated the variational parame- 
ter A using the normalization constraint: 



d.v 



\A\ 2 TT^ 2 W) 



and by requiring that the separate parts of the product 
wave function to be separately normalized to unity: 



dy\4>\ 2 = l, \A\ 2 n 1/2 w = I. (16) 



The second step in the HLVM is to apply the Euler- 
Lagrangc equations of motion to Lhybrid to obtain the 
equations of motion. The equation for <fi is a modified 
version of Eq. (|11|) : 



/ 9Lhybrid\ 


d 


( <9£hybrid\ 


V m ) 


+ 0~y 





dL 



hybrid 



dej)* 



(17) 



and the Euler-Lagrange equations for w and /3 are the 
usual ones: 



d ( 9-Lhybrid 

di \ d~q 



dL 



hybrid 



dq 



= 0, q = w,/3. 



(18) 
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Applying Eq. (fT7|) yields the following equation for 



.d(j> 



_ ( d 2 4> d 2 4> 
di \ dx 2 dy 2 

+ F{T)cp 

where 



V(x,y)(j) 



gN 



2irw 



(19) 



F(i) = \fiw 2 + 2^ + 2P 2 w 2 + \X 2 w 2 . (20) 

This seemingly complicated function of t can be trans- 
formed away by defining 



(21) 



Inserting this into the equation of motion for <f> yields an 
effective 2D Gross-Pitaevskii-likc equation for (ft: 



dej) 
'fit 



d 2 4> d 2 4> 
dx 2 dy 2 



V{x,y)4>- 



gN 



2irw 



2 _ 



(22) 



Applying the Euler-Lagrange equation for /3 gives the 
following result: 



^_|0| 2 ) (i^ 2 ) + |0| 2 - 4^ 2 ) 



We can obtain a simplified equation of motion by inte- 
grating both sides of the above over all x and y: 



dy 



dx 



d_ 

di 



A/3w 2 ) 



dy I dx |0| 

(24) 



It is easy to show that the integral on the left is zero 
by using the equation of motion for (f>. The integral on 
the right is unity by normalization and so we obtain the 
following relationship between j3 and w, w: 



w 
Aw 



(25) 



Thus, if w and w are known, j3 is determined. 

Applying the Euler-Lagrange equation for w gives 



H 2 - £ + A~p 2 w + X 2 wj = \gN 



1 



(26) 

Integrating this equation over all (x,y) on both sides as 
before we obtain 



gNU\\ 



where 



pw + Apw -X + \ 2 w 



U\\$)= / dy / dx\0(x,y,i)\ 



(27) 



(28) 



Note that we have used Eq. (|21j) to replace cj> with 0. 

It is possible to eliminate j3 from the above equation 
by differentiating both sides of Eq. (f25|) with respect to 
time. We obtain 



tw = /3w + (3w = /3w + A(3 2 w, 



(29) 



where the second equality results from using Eq. (|25p to 
replace w with A(3w. Now we see that the right-hand- 
side of the above equation is identical to the first two 
terms on the left-hand-side of Eq. (f27|) . Thus we can 
rewrite this equation as follows: 



iv + AX w = + 



/2/wgNUu 



This is the final equation of motion for w. 



(30) 



C. The HLVM equations of motion and the 
variational initial state 



The full set of HLVM equations of motion consist of a 
2D effective GP-like equation for <f>: 



dej) 



d 2 4> d 2 j> 
dx 2 dy 2 



V(x,y)4>- 



gN 



2irw 



and an equation for w: 



AX w = 



^gNU {l 



2 _ 



(31) 



(32) 



These two equations form a closed system from which 
w(t), w(i), and <j>(x,y,t) can be obtained. From these, 
the value of /3(f) and F(t) can be calculated: 



+ 2(3 2 w 2 + \X 2 w 2 . (33) 



Using these quantities, the full value of the variational 
trial wave function can be constructed: 



Atrial (a, 2/, Z,*) = 



1 



rV2, 



1/2 



x e 



(34) 



where 



dy / dx\(j)(x,y,t)\ A 



(35) 



Note that Eqs. (|3T|) and ([3^)1 are coupled. The nonlincar 
term in the 2D GPE for <f> contains the gaussian width, 
w, while the equation for w contains the factor U\\ which 
is the integral of the fourth power of <f>. 

The final element required for this method to be used 
as a means to find an approximation to the solution of 
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FIG. 1: (color online) A comparison of the vertical column density of a Bose-Einstein condensate, as determined by the 3D 
TDGPE and the HLVM equations of motion, after direct release from a ring-shaped trap is plotted along a line through the 
trap center (x axis) for various times of flight (TOF) during expansion, (a) TOF = 0.0 ms, (b) TOF = 2.0 ms, (c) TOF = 4.0 
ms, (d) TOF = 6.0 ms, (e) TOF = 8.0 ms, and (f) TOF = 10.0 ms. The condensate is given one unit of angular momentum 
before release by phase imprint. 



the 3D TDGPE is a set of initial conditions. We present 
one possibility here based on the physics of Bose-Einstein 
condensate systems. 

In a typical BEC experiment a condensate is formed in 
an atom trap. If no further changes in the condensate's 
environment occur, the condensate wave function should 
then, in principle, only acquire an overall time-dependent 
phase as it evolves in time. In the HLVM this situation 
should therefore be represented by the stationary solu- 
tion of the above equations. We denote this stationary 
solution as 

4>(x, y, 0) = 4>o(x,y) and w{0) = w o . (36) 

This solution satifies the following time-independent 
equations. 



d 2 ^ d 2 o \ - - 
+ I + V{x,y)(/)o 



dx 2 



dy 2 



\y/2lTWoJ 



and 



R(w ) = 4X 2 w 




0. 



where 



^,0 = 



dy / dx\(f> (x,y)\ 4 



(37) 



(38) 



(39) 



The factor /i in the equation for <p is the chemical po- 
tential of the initial condensate. These equations for 
the stationary variational solution must be solved self- 
consistently. 



III. COMPARISON WITH 3D GPE 

In this section we illustrate the ability of the HLVM 
equations of motion to approximate the exact solution of 
the 3D TDGPE by comparing the two solutions for a case 
of current interest. The system we will consider is that 
of a Bose-Einstein condensate of 23 Na atoms confined in 
a ring-shaped potential under the same conditions as in 
a recent experiment conducted at NIST. To simplify 
the analysis we will only compute the profile of the con- 
densate density integrated along the vertical direction for 
each point in the plane of weak confinement. This quan- 
tity predicted by the TDGPE will be compared with that 
predicted by the HLVM equations of motion. This is the 
quantity that can be compared with experiment. 

In the NIST experiment a vertical Laguerre-Gauss 
(LG) laser beam (LGj) was intersected with a horizon- 
tal light sheet. The shape of the vertical LG beam was 
approximately a hollow cylinder with thick walls so that 
its intersection with the horizontal light sheet created 
a ring-shaped region of maximum light intensity. Tun- 
ing the frequency of the beams to the red of the lowest 
electronic transition created a potential that caused the 
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FIG. 2: (color online) GPE/LVM vertical column density comparison for different angular momenta applied to the initial state. 
The condensates are formed, stirred to add m units of angular momentum (simulated by phase imprint), and then released and 
allowed to expand for 10 ms. (a) m = 0, (b) m = 1, (c) m = 2, (d) m = 3, (e) m = 4, and (f) m = 5. 



atoms to seek the maximum intensity. 

In this comparison, wc simulate an experiment in 
which a condensate is created in this ring potential, op- 
tionally stirred, and then probed. We simulate two types 
of probes: (1) direct release of the condensate by turning 
off all trapping potentials after stirring, and (2) release 
of the condensate after the ramp down of the Laguerre- 
Gauss potential. The stirring, which adds m units of 
angular momentum to the condensate, is simulated by 
phase imprint. That is, the initial condensate wave func- 
tion is multiplied by e lm ^ where 4> is the azimuthal angle 
around the vertical z axis and m is an integer. 

In each case we will compare what would be the mea- 
sured density profile, as predicted by the 3D GPE and 
by the HLVM, for different times during the ramp down 
or expansion where the value of m is fixed, and also for a 
fixed final timc-of-flight for a range of different m values. 
For maximum clarity, we present the two density profiles 
as a plot of the density along a line that cuts through 
the center of the ring. Since all of the density profiles are 
cylindrically symmetric, these plots will convey all of the 
available density information. 

In these simulations the trap potential is modeled as 
the sum of a Laguerre-Gauss optical potential [l8| plus 
a vertical gaussian due to the light sheet. This potential 
can be written (in scaled units) as: 

V tiap (x,y,z,t) = -etW(i) (t^) 

- V sheet +\ 2 z 2 = U l{ (x,y) + X 2 z 2 , 

(40) 



where the factor e = 2.718 .. . is included so that Vlg 
becomes the depth of the potential due to the Lagucrre- 
Gaussian beam and Tm is the radial position of its min- 
imum. A time-dependent, dimensionless turn-on func- 
tion, < /(t) < 1, is inserted to simulate the ramp 
down of the Laguerre-Gauss potential. The factor V s h e et 
is the depth of the potential due to the light sheet and 
the z dependent gaussian light sheet potential has been 
approximated by an harmonic oscillator. 

Both the 3D TDGPE and the 2D GPE part of the 
HLVM equations of motion were solved using the split- 
step, Crank-Nicolson method. The 3D TDGPE was 
solved on a grid in which there were 400 points along 
x and y and 200 points along z. The 2D GPE part of the 
HLVM equations of motion was solved on a grid of 800 
points along both x and y. The codes that were used to 
solve these equations were extensively modified versions 
of codes publicly available in the literature [lj| . The ini- 
tial condensate wave function for the 3D TDGPE was 
obtained by solving it in imaginary time. Initial condi- 
tions for the HLVM equation of motion were obtained by 
solving equations ([37]) and (|38|) self consistently as fol- 
lows. First, a value for wq was chosen, the associated <f> 
was then found by integrating Eq. pip in imaginary time, 
next the value of J/y was calculated which was then used 
to compute the value of R(wo) in Eq. P5|) . The value 
of wo was incremented and the process was repeated to 
compute a new value of R. This process was continued 
until a root of R(wq) was found. The value of wq for this 
case is the self-consistent gaussian width of the station- 
ary solution of the HLVM equations of motion. 
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FIG. 3: (color online) GPE/LVM vertical column density comparison for the case where the LG potential ramped down to 20% 
of its initial depth over a span of 50 ms. A phase imprint is applied that simulates one unit of angular momentum added by 
stirring. The plots show the comparison for ramp times of (a) ms, (b) 10 ms, (c) 20 ms, (d) 30 ms, (e) 40 ms, and (f) 50 ms. 



In the direct-release process simulated, the number of 
condensate atoms was N =750,000 atoms and the scat- 
tering length of 23 Na was taken to be 53 bohr. The min- 
imum of the LG potential was set at tm = 24^m. The 
depth of the LG potential was taken to be Vlg = 227 
nK which is equivalent (via Vlg = \Mu>^r\ [ ) to a ra- 
dial harmonic frequency of cj r /(27r) = 120 Hz. The fre- 
quency of the harmonic light sheet potential was taken 
as uj z /(2tt) = 320 Hz and the light-sheet depth was 
Vsheet = 473 nK although this last quantity makes no 
difference in the shape of the initial-state density. 

Figure Q] displays a comparison of the integrated col- 
umn density of a released ring BEC predicted by the 3D 
TDGPE with that predicted by the HLVM equations of 
motion at six different times-of-flight (TOF) after release 
beginning with Fig. (jTJt) showing the moment of release. 
The condensate has been stirred so that it is released hav- 
ing one unit of angular momentum. We note that there 
is good agreement with quantitative differences occuring 
in the heights of individual peaks and in the position of 
the peaks at later times. The comparisons are typical of 
a variational solution in that they are the "best fit" to 
the exact solution for the given trial wave function. 

A comparison of the GPE and LVM results for a ring 
BEC directly released from the trap and allowed to ex- 
pand for a fixed TOF for different initial angular mo- 
menta is exhibited in Fig. [5J The figure displays com- 
parisons for m values ranging from to 5. Again the 
agreement is good although there are some quantitative 
differences as to positions of the individual peaks. It 
is clear that there is qualitative and almost-quantitative 



agreement between the 3D TDGPE and the HLVM equa- 
tions of motion for these cases. 

We next compare the results of the TDGPE and HLVM 
for ring-BEC evolution while the LG potential is ramped 
down from its initial value. This differs from the pre- 
vious comparison in that the confining light -sheet po- 
tential remains unchanged during the ramp down. In 
the simulated ramp-down process the number of con- 
densate atoms was N =500,000 atoms. The LG potential 
depth was ramped linearly down from its initial value of 
Vlg = 227 nK (the same as previously) to 20% of this 
value over a timespan of 50 ms. The light sheet poten- 
tial was the same as in the direct-release simulations. A 
phase imprint was applied to the condensate to simulate 
one unit of angular momentum added by stirring. 

Figure [3] shows the comparison starting at t = ms 
and for every 10 ms thereafter until the rampdown is 
complete. The two solutions again exhibit the type of 
agreement that is usual for variational approximations 
in that the variational solution is a "best fit" to the nu- 
merical solultion. With that caveat, the agreement here 
is quite good during the entire ramp down. 



IV. SUMMARY 

In this paper we derived equations of motion whose 
solution approximates the solution of the 3D TDGPE 
applied to a quasi-2D Bose-Einstein condensate. The 
equations were derived using a hybrid Lagrangian Vari- 
ational Method. Similar equations of motion were de- 
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rived earlier [l7| for quasi-lD BEC systems. The main 
advantage of solving these equations is that numerical 
solution of the HLVM equation can be performed 100 to 
1000 times faster than solving the 3D TDGPE. In the 
comparison simulation presented in Section IIII1 solving 
the realtime 3D TDGPE required more than 24 hours of 
CPU time while solving the HLVM equation took about 
10 minutes on a commodity desktop PC. The resulting 
speedup here is roughly a factor of 150. 

This advantage enables rapid simulation of many dif- 
ferent possible quasi-2D systems. It should be noted 
that when the HLVM was applied to a quasi-lD sys- 
tem describing soliton splitting in optical fibers [2(j, it 
was found that occasionally the HLVM equations of mo- 
tion did not provide any advantage over the regular LVM 
technique. In the work cited, the authors recommended 
that any important results that come out of HLVM sim- 
ulations be confirmed by simulations using the full equa- 
tions. We did not find any case where the HLVM equa- 
tions predicted behavior that was qualitatively different 
from the 3D TDGPE. However, we agree with the rec- 



ommendation of the authors of Ref. [2C§ . 

This caveat notwithstanding, the HLVM equations de- 
rived in this paper enable rapid study of different sys- 
tems of current experimental interest. In particular, they 
should be useful in simulating time-dependent behavior 
of quasi-2D atomtronic systems where mean-field the- 
ory applies. We expect this approximation to become 
a useful tool in studying future quasi-2D Bosc-Einstein 
condensate systems. 
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